library(xtable)
library(stargazer)
library(multiwayvcov)
library(lmtest)
library(scales)
library(foreign)
library(nFactors)
library(FactoMineR)
rm(list = ls())

#####upload data#####
load(file='~/Dropbox/many ways to be right submission/BJPS final submission/BJPS replication materials/data_2017.RData')
load(file='~/Dropbox/many ways to be right submission/BJPS final submission/BJPS replication materials/data_1990.RData')

#####Table 1######
share.cross.1990 <- 100*(table(data.1990$quadrants) / nrow(data.1990))
share.cross.2017 <- 100*(table(data.2017$quadrants) / nrow(data.2017))
share.cross.1990 <- round(share.cross.1990, digits = 0)
share.cross.2017 <- round(share.cross.2017, digits = 0)
share.cross <- rbind(share.cross.1990, share.cross.2017)
rownames(share.cross) <- c("1990 sample", "2017 sample")
colnames(share.cross) <- c("Consistent progressives", "Consistent conservatives", "Market cosmopolitans", "welfare chauvinists")
xtable(share.cross)

######Table 2########

right_cont.quadrants.90 <- lm(lr_cont~quadrants + as.factor(countries), 
                              weight=weights, 
                              data=data.1990)
right_cont.quadrants.90.clustered <- coeftest(right_cont.quadrants.90, 
                                              cluster.vcov(right_cont.quadrants.90, data.1990$countries))

right_cont.quadrants.cov.90 <- lm(lr_cont~quadrants + gender + age + as.factor(income) + rel.or + 
                                    unions + as.factor(countries), 
                                  weight=weights, 
                                  data=data.1990)
right_cont.quadrants.cov.90.clustered <- coeftest(right_cont.quadrants.cov.90, 
                                                  cluster.vcov(right_cont.quadrants.cov.90, data.1990$countries))

right_b.quadrants.90 <- lm(binary.right~quadrants + as.factor(countries), 
                           weight=weights, 
                           data=data.1990)
right_b.quadrants.90.clustered <- coeftest(right_b.quadrants.90, 
                                           cluster.vcov(right_b.quadrants.90, data.1990$countries))

right_b.quadrants.cov.90 <- lm(binary.right~quadrants + gender + age + as.factor(income) + rel.or + 
                                 unions + as.factor(countries), 
                               weight=weights, 
                               data=data.1990)
right_b.quadrants.cov.90.clustered <- coeftest(right_b.quadrants.cov.90, 
                                               cluster.vcov(right_b.quadrants.cov.90, data.1990$countries))

right_cont.quadrants.17 <- lm(lr_cont~quadrants + as.factor(countries), 
                              weight=weights, 
                              data=data.2017)
right_cont.quadrants.17.clustered <- coeftest(right_cont.quadrants.17, 
                                              cluster.vcov(right_cont.quadrants.17, data.2017$countries))

right_cont.quadrants.cov.17 <- lm(lr_cont~quadrants + gender + age + as.factor(income) + high.ed + rel.or + 
                                    unions  + as.factor(countries), 
                                  weight=weights, 
                                  data=data.2017)
right_cont.quadrants.cov.17.clustered <- coeftest(right_cont.quadrants.cov.17, 
                                                  cluster.vcov(right_cont.quadrants.cov.17, data.2017$countries))

right_b.quadrants.17 <- lm(binary.right~quadrants + as.factor(countries), 
                           weight=weights, 
                           data=data.2017)
right_b.quadrants.17.clustered <- coeftest(right_b.quadrants.17, 
                                           cluster.vcov(right_b.quadrants.17, data.2017$countries))

right_b.quadrants.cov.17 <- lm(binary.right~quadrants + gender + age + as.factor(income) + high.ed + rel.or + 
                                 unions  + as.factor(countries), 
                               weight=weights, 
                               data=data.2017)
right_b.quadrants.cov.17.clustered <- coeftest(right_b.quadrants.cov.17, 
                                               cluster.vcov(right_b.quadrants.cov.17, data.2017$countries))

stargazer(right_cont.quadrants.90.clustered, right_cont.quadrants.cov.90.clustered,
          right_cont.quadrants.17.clustered, right_cont.quadrants.cov.17.clustered,
          right_b.quadrants.90.clustered, right_b.quadrants.cov.90.clustered,
          right_b.quadrants.17.clustered, right_b.quadrants.cov.17.clustered,
          omit="countries", no.space=TRUE,
          covariate.labels = c("Consistent conservatives", "Market cosmopolitans", "Welfare chauvinists",
                               "Gender (1=female)", "Age", "Income: medium", "Income: high",
                               "Higher education (=1)",  "Member in religious organization (=1)", "Union member (=1)"))

#for number of observations in R^2
stargazer(right_cont.quadrants.90, right_cont.quadrants.cov.90,
          right_cont.quadrants.17, right_cont.quadrants.cov.17,
          right_b.quadrants.90, right_b.quadrants.cov.90,
          right_b.quadrants.17, right_b.quadrants.cov.17)

######Table 3########

right_cont.delta.90 <- lm(lr_cont~red.minus.values + I(red.minus.values^2) + as.factor(countries), 
                          weight=weights, 
                          data=data.1990)
right_cont.delta.90.clustered <- coeftest(right_cont.delta.90, 
                                          cluster.vcov(right_cont.delta.90, data.1990$countries))

right_cont.delta.covariates.90 <- lm(lr_cont~red.minus.values + I(red.minus.values^2) + 
                                       gender + age + as.factor(income) + rel.or + unions  + as.factor(countries), 
                                     weight=weights, 
                                     data=data.1990)
right_cont.delta.covariates.90.clustered <- coeftest(right_cont.delta.covariates.90, 
                                                     cluster.vcov(right_cont.delta.covariates.90, data.1990$countries))

right_b.delta.90 <- lm(binary.right~red.minus.values + I(red.minus.values^2) + 
                         as.factor(countries), 
                       weight=weights,
                       data=data.1990)
right_b.delta.90.clustered <- coeftest(right_b.delta.90, 
                                       cluster.vcov(right_b.delta.90, data.1990$countries))

right_b.delta.covariates.90 <- lm(binary.right~red.minus.values + I(red.minus.values^2) + 
                                    gender + age + as.factor(income) + rel.or + unions + as.factor(countries), 
                                  weight=weights,
                                  data=data.1990)
right_b.delta.covariates.90.clustered <- coeftest(right_b.delta.covariates.90, 
                                                  cluster.vcov(right_b.delta.covariates.90, data.1990$countries))

right_cont.delta.17 <- lm(lr_cont~red.minus.values + I(red.minus.values^2) + as.factor(countries), 
                          weight=weights, 
                          data=data.2017)
right_cont.delta.17.clustered <- coeftest(right_cont.delta.17, 
                                          cluster.vcov(right_cont.delta.17, data.2017$countries))

right_cont.delta.covariates.17 <- lm(lr_cont~red.minus.values + I(red.minus.values^2) + 
                                       gender + age + as.factor(income) + high.ed + rel.or + unions  + as.factor(countries), 
                                     weight=weights, 
                                     data=data.2017)
right_cont.delta.covariates.17.clustered <- coeftest(right_cont.delta.covariates.17, 
                                                     cluster.vcov(right_cont.delta.covariates.17, data.2017$countries))

right_b.delta.17 <- lm(binary.right~red.minus.values + I(red.minus.values^2) + as.factor(countries), 
                       weight=weights, 
                       data=data.2017)
right_b.delta.17.clustered <- coeftest(right_b.delta.17, 
                                       cluster.vcov(right_b.delta.17, data.2017$countries))

right_b.delta.covariates.17 <- lm(binary.right~red.minus.values + I(red.minus.values^2) + 
                                    gender + age + as.factor(income) + high.ed + rel.or + unions  + as.factor(countries), 
                                  weight=weights, 
                                  data=data.2017)
right_b.delta.covariates.17.clustered <- coeftest(right_b.delta.covariates.17, 
                                                  cluster.vcov(right_b.delta.covariates.17, data.2017$countries))

stargazer(right_cont.delta.90.clustered, 
          right_cont.delta.covariates.90.clustered,
          right_cont.delta.17.clustered,
          right_cont.delta.covariates.17.clustered,
          right_b.delta.90.clustered, 
          right_b.delta.covariates.90.clustered,
          right_b.delta.17.clustered, 
          right_b.delta.covariates.17.clustered,
          omit="countries", no.space=TRUE,
          covariate.labels = c("Delta EC", "Delta EC^2", "Gender (1=female)", "Age", "Income: medium", "Income: High",
                               "Higher education (=1)",  "Member in religious organization (=1)", "Union member (=1)"))

stargazer(right_cont.delta.90, 
          right_cont.delta.covariates.90,
          right_cont.delta.17,
          right_cont.delta.covariates.17,
          right_b.delta.90, 
          right_b.delta.covariates.90,
          right_b.delta.17, 
          right_b.delta.covariates.17)

#######Figure 2########

delta.x.90 <- seq(min(as.numeric(data.1990$red.minus.values), na.rm=T),
                  max(as.numeric(data.1990$red.minus.values), na.rm=T),0.1)
xnew.delta.90=list(red.minus.values=delta.x.90, 
                   gender=rep(0, length(delta.x.90)), 
                   age=rep(42, length(delta.x.90)),
                   income=rep(2, length(delta.x.90)),
                   rel.or=rep(0, length(delta.x.90)),
                   unions=rep(0, length(delta.x.90)),
                   countries=rep("Great Britain",length(delta.x.90)),
                   length(delta.x.90))
xnew.delta.90.nocov=list(red.minus.values=delta.x.90, 
                         countries=rep("Great Britain",length(delta.x.90)),
                         length(delta.x.90))

pred_b.delta.90 <-predict(right_b.delta.covariates.90, newdata=xnew.delta.90, 
                          type="response", se=TRUE)
pred_cont.delta.90 <-predict(right_cont.delta.covariates.90, newdata=xnew.delta.90, 
                             type="response", se=TRUE)

pred_b.delta.nocov.90 <-predict(right_b.delta.90, newdata=xnew.delta.90.nocov, 
                                type="response", se=TRUE)

pred_cont.delta.nocov.90 <-predict(right_cont.delta.90, newdata=xnew.delta.90.nocov, 
                                   type="response", se=TRUE)

delta.x.17 <- seq(min(as.numeric(data.2017$red.minus.values), na.rm=T),
                  max(as.numeric(data.2017$red.minus.values), na.rm=T),0.1)
xnew.delta.17=list(red.minus.values=delta.x.17, 
                   gender=rep(0, length(delta.x.17)), 
                   age=rep(42, length(delta.x.17)),
                   income=rep(2, length(delta.x.17)),
                   rel.or=rep(0, length(delta.x.17)),
                   high.ed=rep(0, length(delta.x.17)),
                   unions=rep(0, length(delta.x.17)),
                   countries=rep("Great Britain",length(delta.x.17)),
                   length(delta.x.17))
xnew.delta.17.nocov=list(red.minus.values=delta.x.17, 
                         countries=rep("Great Britain",length(delta.x.17)),
                         length(delta.x.17))

pred_b.delta.17 <-predict(right_b.delta.covariates.17, newdata=xnew.delta.17, 
                          type="response", se=TRUE)
pred_cont.delta.17 <-predict(right_cont.delta.covariates.17, newdata=xnew.delta.17, 
                             type="response", se=TRUE)

pred_b.delta.nocov.17 <-predict(right_b.delta.17, newdata=xnew.delta.17.nocov, 
                                type="response", se=TRUE)
pred_cont.delta.nocov.17 <-predict(right_cont.delta.17, newdata=xnew.delta.17.nocov, 
                                   type="response", se=TRUE)

#pdf(file='right_b_delta_90_17.pdf') 
par(mar=c(5.1,4.6,4.1,2.1))
plot(delta.x.17, pred_cont.delta.17$fit, type="n", ylim=c(0, 01), xlim=c(-6, 6),
     cex.axis=1.3, cex.lab=1.5, ylab="Identifying with the Right=1",
     xlab=expression(paste(Delta, "EC")))

axis(1, at=c(-70, 70), lab=c("Welfare \n chauvinists", 
                             "Market \n cosmopolitans"),
     cex.axis=1.5, tick = FALSE, line=3)

lines(delta.x.17, pred_b.delta.17$fit, lwd=6, col="blue", lty=5)
lines(delta.x.17, pred_b.delta.17$fit, lwd=6, col="blue", lty=5)

lines(delta.x.90, pred_b.delta.90$fit, lwd=6, col="red", lty=1)

legend("topright", c("1990", "2017"), cex=2, col=c("red", "blue"), lty=c(1,5), lwd=c(6,6), bty="n")
rug(data.1990$red.minus.values, col="red")
rug(data.2017$red.minus.values, side=3, col="blue")
#dev.off()

#####Table A2########

lm.economic.1990 <- lm(econ_pref~gender + age + as.factor(income) + rel.or + unions  + as.factor(countries), 
                       weight=weights, 
                       data=data.1990)
lm.economic.1990.clustered <- coeftest(lm.economic.1990, 
                                       cluster.vcov(lm.economic.1990, data.1990$countries))

lm.cultural.1990 <- lm(cultural_pref~gender + age + as.factor(income) + rel.or + unions  + as.factor(countries), 
                       weight=weights, 
                       data=data.1990)
lm.cultural.1990.clustered <- coeftest(lm.cultural.1990, 
                                       cluster.vcov(lm.cultural.1990, data.1990$countries))

lm.economic.2017 <- lm(econ_pref~gender + age + as.factor(income) + high.ed + rel.or + unions + as.factor(countries), 
                       weight=weights, 
                       data=data.2017)
lm.economic.2017.clustered <- coeftest(lm.economic.2017, 
                                       cluster.vcov(lm.economic.2017, data.2017$countries))

lm.cultural.2017 <- lm(cultural_pref~gender + age + as.factor(income) + high.ed + rel.or + unions  + as.factor(countries), 
                       weight=weights, 
                       data=data.2017)
lm.cultural.2017.clustered <- coeftest(lm.cultural.2017, 
                                       cluster.vcov(lm.cultural.2017, data.2017$countries))

stargazer(lm.economic.1990.clustered, lm.cultural.1990.clustered,
          lm.economic.2017.clustered, lm.cultural.2017.clustered,
          omit="countries", no.space=TRUE,
          covariate.labels = c("Gender (1=female)", "Age", "Income (medium)" , "Income (high)",
                               "Highed education (=1)", "Member in religious organization (=1)",  
                               "Union member (=1)"))

stargazer(lm.economic.1990, lm.cultural.1990,
          lm.economic.2017, lm.cultural.2017,
          omit="countries", no.space=TRUE,
          covariate.labels = c("Gender (1=female)", "Age", "Income (medium)" , "Income (high)",
                               "Highed education (=1)", "Member in religious organization (=1)",  
                               "Union member (=1)"))

####Table A3####

welfare.chauvinists.lm <- lm(welfare_chauvinists_voters~ gender + age + as.factor(income) + high.ed + rel.or + unions  + as.factor(countries), 
                             weight=weights, 
                             data=data.2017)
welfare.chauvinists.lm.clustered <- coeftest(welfare.chauvinists.lm, cluster.vcov(welfare.chauvinists.lm, data.2017$countries))

market.cosmopolitans.lm <- lm(market_cosmopolitans_voters~ gender + age + as.factor(income) + high.ed + rel.or + unions  + as.factor(countries), 
                              weight=weights, 
                              data=data.2017)
market.cosmopolitans.lm.clustered <- coeftest(market.cosmopolitans.lm, cluster.vcov(market.cosmopolitans.lm, data.2017$countries))

stargazer(welfare.chauvinists.lm.clustered,market.cosmopolitans.lm.clustered,
          omit="countries", no.space=TRUE,
          covariate.labels = c("Gender (1=female)", "Age", "Income: medium", "Income: high",
                               "Higher education (=1)",  "Member in religious organization (=1)", "Union member (=1)"))
stargazer(welfare.chauvinists.lm ,market.cosmopolitans.lm)

#######Table A4######

work.welfare <- lm(work.important~ welfare_chauvinists_voters + gender + age + high.ed + as.factor(income) + rel.or + 
                     unions  + as.factor(countries), 
                   weight=weights, 
                   data=data.2017)
work.welfare.clustered <- coeftest(work.welfare, cluster.vcov(work.welfare, data.2017$countries))

work.market <- lm(work.important~ market_cosmopolitans_voters + gender + age + high.ed + as.factor(income) + rel.or + 
                    unions  + as.factor(countries), 
                  weight=weights, 
                  data=data.2017)
work.market.clustered <- coeftest(work.market, cluster.vcov(work.market, data.2017$countries))

religion.welfare <- lm(religion.important~ welfare_chauvinists_voters + gender + age + high.ed + as.factor(income) + rel.or + 
                         unions  + as.factor(countries), 
                       weight=weights, 
                       data=data.2017)
religion.welfare.clustered <- coeftest(religion.welfare, cluster.vcov(religion.welfare, data.2017$countries))

religion.market <- lm(religion.important~ market_cosmopolitans_voters + gender + age + high.ed + as.factor(income) + rel.or + 
                        unions  + as.factor(countries), 
                      weight=weights, 
                      data=data.2017)
religion.market.clustered <- coeftest(religion.market, cluster.vcov(religion.market, data.2017$countries))

stargazer(work.welfare.clustered, work.market.clustered,
          religion.welfare.clustered, religion.market.clustered,
          omit="countries", no.space=TRUE,
          covariate.labels = c("Welfare chauvinists", "Market cosmopolitans",
                               "Gender (1=female)", "Age", "Higher education", "Income: medium", "Income: high",
                               "Member in religious organization (=1)", "Union member (=1)"))
stargazer(work.welfare, work.market,
          religion.welfare, religion.market)

####Tables A5-A11######

countries.christian.democrats <- data.2017[data.2017$countries %in% c("Austria", "Finland", "Germany", "Italy", "Netherlands", "Norway"),]
countries.conservatives <- data.2017[data.2017$countries %in% c("Finland", "France","Iceland", "Italy", "Norway", 
                                                                                    "Spain", "Great Britain"),]
countries.liberals <- data.2017[data.2017$countries %in% c("Austria", "France","Germany", "Iceland", "Italy", "Netherlands", "Norway", 
                                                                               "Spain", "Great Britain"),]
countries.radical.right <- data.2017[data.2017$countries %in% c("Austria", "Finland", "France","Germany", "Italy", "Netherlands", "Norway", 
                                                                                    "Great Britain"),]

countries.social.democrats <- data.2017[data.2017$countries %in% c("Austria", "Finland", "France","Germany", "Iceland", "Italy", "Netherlands", "Norway", 
                                                                                       "Spain", "Great Britain"),]
countries.greens <- data.2017[data.2017$countries %in% c("Austria", "Finland", "France","Germany", "Iceland", "Netherlands", "Norway", "Great Britain"),]
countries.radical.left <- data.2017[data.2017$countries %in% c("Finland", "France", "Germany", "Italy", "Netherlands", "Norway", "Spain"),]

lm.radical_right.welfare <- lm(radical_right~ welfare_chauvinists_voters + 
                                 gender + age + as.factor(income) + high.ed + rel.or + unions  + as.factor(countries), 
                               weight=weights, 
                               data=countries.radical.right)
lm.radical_right.welfare.clustered <- coeftest(lm.radical_right.welfare, 
                                               cluster.vcov(lm.radical_right.welfare, countries.radical.right$countries))

lm.radical_right.market <- lm(radical_right~ market_cosmopolitans_voters +
                                gender + age + as.factor(income) + high.ed + rel.or + unions  + as.factor(countries), 
                              weight=weights, 
                              data=countries.radical.right)
lm.radical_right.market.clustered <- coeftest(lm.radical_right.market, 
                                              cluster.vcov(lm.radical_right.market, countries.radical.right$countries))

lm.radical_right.conservative <- lm(radical_right~ conservative_voters +
                                      gender + age + as.factor(income) + high.ed + rel.or + unions  + as.factor(countries), 
                                    weight=weights, 
                                    data=countries.radical.right)
lm.radical_right.conservative.clustered <- coeftest(lm.radical_right.conservative, 
                                                    cluster.vcov(lm.radical_right.conservative, countries.radical.right$countries))

lm.radical_right.progressive <- lm(radical_right~ progressive_voters +
                                     gender + age + as.factor(income) + high.ed + rel.or + unions  + as.factor(countries), 
                                   weight=weights, 
                                   data=countries.radical.right)
lm.radical_right.progressive.clustered <- coeftest(lm.radical_right.progressive, 
                                                   cluster.vcov(lm.radical_right.progressive, countries.radical.right$countries))

stargazer(lm.radical_right.welfare.clustered,lm.radical_right.market.clustered,
          lm.radical_right.conservative.clustered,lm.radical_right.progressive.clustered,
          omit="countries", no.space=TRUE,
          covariate.labels = c("Welfare chauvinists", "Market cosmopolitans",  "Consistent conservatives", "Consistent progressives",
                               "Gender (1=female)", "Age", "Income: medium", "Income: high",
                               "Higher education (=1)",  "Member in religious organization (=1)", "Union member (=1)"))
stargazer(lm.radical_right.welfare, lm.radical_right.market,
          lm.radical_right.conservative, lm.radical_right.progressive)

rr.coef <- c(lm.radical_right.progressive.clustered[2,1], lm.radical_right.conservative.clustered[2,1],
             lm.radical_right.market.clustered[2,1],lm.radical_right.welfare.clustered[2,1])

rr.sd <- c(lm.radical_right.progressive.clustered[2,2], lm.radical_right.conservative.clustered[2,2],
           lm.radical_right.market.clustered[2,2],lm.radical_right.welfare.clustered[2,2])

##

lm.christian_democrats.welfare_chuavinists <- lm(christian_democrats~ welfare_chauvinists_voters + 
                                                   gender + age + as.factor(income) + high.ed + rel.or + unions  + as.factor(countries), 
                                                 weight=weights, 
                                                 data=countries.christian.democrats)
lm.christian_democrats.welfare_chuavinists.clustered <- coeftest(lm.christian_democrats.welfare_chuavinists, 
                                                                 cluster.vcov(lm.christian_democrats.welfare_chuavinists, countries.christian.democrats$countries))

lm.christian_democrats.conservative_voters <- lm(christian_democrats~ conservative_voters +  
                                                   gender + age + as.factor(income) + high.ed + rel.or + unions  + as.factor(countries), 
                                                 weight=weights, 
                                                 data=countries.christian.democrats)
lm.christian_democrats.conservative_voters.clustered <- coeftest(lm.christian_democrats.conservative_voters, 
                                                                 cluster.vcov(lm.christian_democrats.conservative_voters, countries.christian.democrats$countries))

lm.christian_democrats.market_cosmopolitans_voters <- lm(christian_democrats~ market_cosmopolitans_voters +  
                                                           gender + age + as.factor(income) + high.ed + rel.or + unions  + as.factor(countries), 
                                                         weight=weights, 
                                                         data=countries.christian.democrats)
lm.christian_democrats.market_cosmopolitans_voters.clustered <- coeftest(lm.christian_democrats.market_cosmopolitans_voters, 
                                                                         cluster.vcov(lm.christian_democrats.market_cosmopolitans_voters, countries.christian.democrats$countries))

lm.christian_democrats.progressives <- lm(christian_democrats~ progressive_voters +  
                                            gender + age + as.factor(income) + high.ed + rel.or + unions  + as.factor(countries), 
                                          weight=weights, 
                                          data=countries.christian.democrats)
lm.christian_democrats.progressives.clustered <- coeftest(lm.christian_democrats.progressives, 
                                                          cluster.vcov(lm.christian_democrats.progressives, countries.christian.democrats$countries))

stargazer(lm.christian_democrats.welfare_chuavinists.clustered, lm.christian_democrats.market_cosmopolitans_voters.clustered,
          lm.christian_democrats.conservative_voters.clustered, lm.christian_democrats.progressives.clustered,
          omit="countries", no.space=TRUE,
          covariate.labels = c("Welfare chauvinists", "Market cosmopolitans",  "Consistent conservatives", "Consistent progressives",
                               "Gender (1=female)", "Age", "Income: medium", "Income: high",
                               "Higher education (=1)",  "Member in religious organization (=1)", "Union member (=1)"))

stargazer(lm.christian_democrats.welfare_chuavinists, lm.christian_democrats.market_cosmopolitans_voters,
          lm.christian_democrats.conservative_voters, lm.christian_democrats.progressives)

cd.coef <- c(lm.christian_democrats.progressives.clustered[2,1], lm.christian_democrats.conservative_voters.clustered[2,1],
             lm.christian_democrats.market_cosmopolitans_voters.clustered[2,1],lm.christian_democrats.welfare_chuavinists.clustered[2,1])

cd.sd <- c(lm.christian_democrats.progressives.clustered[2,2], lm.christian_democrats.conservative_voters.clustered[2,2],
           lm.christian_democrats.market_cosmopolitans_voters.clustered[2,2],lm.christian_democrats.welfare_chuavinists.clustered[2,2])

#

lm.conservatives.welfare_chauvinists_voters <- lm(conservatives~ welfare_chauvinists_voters +  
                                                    gender + age + as.factor(income) + high.ed + rel.or + unions  + as.factor(countries), 
                                                  weight=weights, 
                                                  data=countries.conservatives)
lm.conservatives.welfare_chauvinists_voters.clustered <- coeftest(lm.conservatives.welfare_chauvinists_voters, 
                                                                  cluster.vcov(lm.conservatives.welfare_chauvinists_voters, countries.conservatives$countries))

lm.conservatives.market_cosmopolitan_voters <- lm(conservatives~ market_cosmopolitans_voters +  
                                                    gender + age + as.factor(income) + high.ed + rel.or + unions  + as.factor(countries), 
                                                  weight=weights, 
                                                  data=countries.conservatives)
lm.conservatives.market_cosmopolitan_voters.clustered <- coeftest(lm.conservatives.market_cosmopolitan_voters, 
                                                                  cluster.vcov(lm.conservatives.market_cosmopolitan_voters, countries.conservatives$countries))

lm.conservatives.conservative_voters <- lm(conservatives~ conservative_voters +  
                                             gender + age + as.factor(income) + high.ed + rel.or + unions  + as.factor(countries), 
                                           weight=weights, 
                                           data=countries.conservatives)
lm.conservatives.conservative_voters.clustered <- coeftest(lm.conservatives.conservative_voters, 
                                                           cluster.vcov(lm.conservatives.conservative_voters, countries.conservatives$countries))

lm.conservatives.progressive_voters <- lm(conservatives~ progressive_voters +  
                                            gender + age + as.factor(income) + high.ed + rel.or + unions  + as.factor(countries), 
                                          weight=weights, 
                                          data=countries.conservatives)
lm.conservatives.progressive_voters.clustered <- coeftest(lm.conservatives.progressive_voters, 
                                                          cluster.vcov(lm.conservatives.progressive_voters, countries.conservatives$countries))

stargazer(lm.conservatives.welfare_chauvinists_voters.clustered, lm.conservatives.market_cosmopolitan_voters.clustered,
          lm.conservatives.conservative_voters.clustered, lm.conservatives.progressive_voters.clustered,
          omit="countries", no.space=TRUE,
          covariate.labels = c("Welfare chauvinists", "Market cosmopolitans",  "Consistent conservatives", "Consistent progressives",
                               "Gender (1=female)", "Age", "Income: medium", "Income: high",
                               "Higher education (=1)",  "Member in religious organization (=1)", "Union member (=1)"))

stargazer(lm.conservatives.welfare_chauvinists_voters, lm.conservatives.market_cosmopolitan_voters,
          lm.conservatives.conservative_voters, lm.conservatives.progressive_voters)

cons.coef <- c(lm.conservatives.progressive_voters.clustered[2,1], lm.conservatives.conservative_voters.clustered[2,1],
               lm.conservatives.market_cosmopolitan_voters.clustered[2,1],lm.conservatives.welfare_chauvinists_voters.clustered[2,1])

cons.sd <- c(lm.conservatives.progressive_voters.clustered[2,2], lm.conservatives.conservative_voters.clustered[2,2],
             lm.conservatives.market_cosmopolitan_voters.clustered[2,2],lm.conservatives.welfare_chauvinists_voters.clustered[2,2])

#
lm.liberals.welfare_cuahvinists_vote <- lm(liberals~  welfare_chauvinists_voters + 
                                             gender + age + as.factor(income) + high.ed + rel.or + unions  + as.factor(countries), 
                                           weight=weights, 
                                           data=countries.liberals)
lm.liberals.welfare_cuahvinists_vote.clustered <- coeftest(lm.liberals.welfare_cuahvinists_vote, 
                                                           cluster.vcov(lm.liberals.welfare_cuahvinists_vote, countries.liberals$countries))

lm.liberals.market_cosmopolitan_vote <- lm(liberals~ market_cosmopolitans_voters +  
                                             gender + age + as.factor(income) + high.ed + rel.or + unions  + as.factor(countries), 
                                           weight=weights, 
                                           data=countries.liberals)
lm.liberals.market_cosmopolitan_vote.clustered <- coeftest(lm.liberals.market_cosmopolitan_vote, 
                                                           cluster.vcov(lm.liberals.market_cosmopolitan_vote, countries.liberals$countries))

lm.liberals.conservative_vote <- lm(liberals~ conservative_voters +  
                                      gender + age + as.factor(income) + high.ed + rel.or + unions  + as.factor(countries), 
                                    weight=weights, 
                                    data=countries.liberals)
lm.liberals.conservative_vote.clustered <- coeftest(lm.liberals.conservative_vote, 
                                                    cluster.vcov(lm.liberals.conservative_vote, countries.liberals$countries))

lm.liberals.progressive_vote <- lm(liberals~ progressive_voters +  
                                     gender + age + as.factor(income) + high.ed + rel.or + unions  + as.factor(countries), 
                                   weight=weights, 
                                   data=countries.liberals)
lm.liberals.progressive_vote.clustered <- coeftest(lm.liberals.progressive_vote, 
                                                   cluster.vcov(lm.liberals.progressive_vote, countries.liberals$countries))

stargazer(lm.liberals.welfare_cuahvinists_vote.clustered, lm.liberals.market_cosmopolitan_vote.clustered,
          lm.liberals.conservative_vote.clustered, lm.liberals.progressive_vote.clustered,
          omit="countries", no.space=TRUE,
          covariate.labels = c("Welfare chauvinists", "Market cosmopolitans",  "Consistent conservatives", "Consistent progressives",
                               "Gender (1=female)", "Age", "Income: medium", "Income: high",
                               "Higher education (=1)",  "Member in religious organization (=1)", "Union member (=1)"))

stargazer(lm.liberals.welfare_cuahvinists_vote, lm.liberals.market_cosmopolitan_vote,
          lm.liberals.conservative_vote, lm.liberals.progressive_vote)

lib.coef <- c(lm.liberals.progressive_vote.clustered[2,1], lm.liberals.conservative_vote.clustered[2,1],
              lm.liberals.market_cosmopolitan_vote.clustered[2,1],lm.liberals.welfare_cuahvinists_vote.clustered[2,1])

lib.sd <- c(lm.liberals.progressive_vote.clustered[2,2], lm.liberals.conservative_vote.clustered[2,2],
            lm.liberals.market_cosmopolitan_vote.clustered[2,2],lm.liberals.welfare_cuahvinists_vote.clustered[2,2])

#

lm.social_democrats.welfare <- lm(social_democrats~welfare_chauvinists_voters +  
                                    gender + age + as.factor(income) + high.ed + rel.or + unions  + as.factor(countries), 
                                  weight=weights, 
                                  data=countries.social.democrats)
lm.social_democrats.welfare.clustered <- coeftest(lm.social_democrats.welfare, 
                                                  cluster.vcov(lm.social_democrats.welfare, countries.social.democrats$countries))

lm.social_democrats.market <- lm(social_democrats~market_cosmopolitans_voters +  
                                   gender + age + as.factor(income) + high.ed + rel.or + unions  + as.factor(countries), 
                                 weight=weights, 
                                 data=countries.social.democrats)
lm.social_democrats.market.clustered <- coeftest(lm.social_democrats.market, 
                                                 cluster.vcov(lm.social_democrats.market, countries.social.democrats$countries))

lm.social_democrats.conservative <- lm(social_democrats~conservative_voters +  
                                         gender + age + as.factor(income) + high.ed + rel.or + unions  + as.factor(countries), 
                                       weight=weights, 
                                       data=countries.social.democrats)
lm.social_democrats.conservative.clustered <- coeftest(lm.social_democrats.conservative, 
                                                       cluster.vcov(lm.social_democrats.conservative, countries.social.democrats$countries))

lm.social_democrats.progressive <- lm(social_democrats~progressive_voters +  
                                        gender + age + as.factor(income) + high.ed + rel.or + unions  + as.factor(countries), 
                                      weight=weights, 
                                      data=countries.social.democrats)
lm.social_democrats.progressive.clustered <- coeftest(lm.social_democrats.progressive, 
                                                      cluster.vcov(lm.social_democrats.progressive, countries.social.democrats$countries))

stargazer(lm.social_democrats.welfare.clustered, lm.social_democrats.market.clustered,
          lm.social_democrats.conservative.clustered, lm.social_democrats.progressive.clustered,
          omit="countries", no.space=TRUE,
          covariate.labels = c("Welfare chauvinists", "Market cosmopolitans",  "Consistent conservatives", "Consistent progressives",
                               "Gender (1=female)", "Age", "Income: medium", "Income: high",
                               "Higher education (=1)",  "Member in religious organization (=1)", "Union member (=1)"))

stargazer(lm.social_democrats.welfare, lm.social_democrats.market,
          lm.social_democrats.conservative, lm.social_democrats.progressive)

socdem.coef <- c(lm.social_democrats.progressive.clustered[2,1], lm.social_democrats.conservative.clustered[2,1],
                 lm.social_democrats.market.clustered[2,1],lm.social_democrats.welfare.clustered[2,1])

socdem.sd <- c(lm.social_democrats.progressive.clustered[2,2], lm.social_democrats.conservative.clustered[2,2],
               lm.social_democrats.market.clustered[2,2],lm.social_democrats.welfare.clustered[2,2])

#
lm.greens.welfare <- lm(greens~welfare_chauvinists_voters +  
                          gender + age + as.factor(income) + high.ed + rel.or + unions  + as.factor(countries), 
                        weight=weights, 
                        data=countries.greens)
lm.greens.welfare.clustered <- coeftest(lm.greens.welfare, 
                                        cluster.vcov(lm.greens.welfare, countries.greens$countries))

lm.greens.market <- lm(greens~market_cosmopolitans_voters +  
                         gender + age + as.factor(income) + high.ed + rel.or + unions  + as.factor(countries), 
                       weight=weights, 
                       data=countries.greens)
lm.greens.market.clustered <- coeftest(lm.greens.market, 
                                       cluster.vcov(lm.greens.market, countries.greens$countries))

lm.greens.conservative <- lm(greens~conservative_voters +  
                               gender + age + as.factor(income) + high.ed + rel.or + unions  + as.factor(countries), 
                             weight=weights, 
                             data=countries.greens)
lm.greens.conservative.clustered <- coeftest(lm.greens.conservative, 
                                             cluster.vcov(lm.greens.conservative, countries.greens$countries))

lm.greens.progressive <- lm(greens~progressive_voters +  
                              gender + age + as.factor(income) + high.ed + rel.or + unions  + as.factor(countries), 
                            weight=weights, 
                            data=countries.greens)
lm.greens.progressive.clustered <- coeftest(lm.greens.progressive, 
                                            cluster.vcov(lm.greens.progressive, countries.greens$countries))

stargazer(lm.greens.welfare.clustered, lm.greens.market.clustered,
          lm.greens.conservative.clustered, lm.greens.progressive.clustered,
          omit="countries", no.space=TRUE,
          covariate.labels = c("Welfare chauvinists", "Market cosmopolitans",  "Consistent conservatives", "Consistent progressives",
                               "Gender (1=female)", "Age", "Income: medium", "Income: high",
                               "Higher education (=1)",  "Member in religious organization (=1)", "Union member (=1)"))

stargazer(lm.greens.welfare, lm.greens.market,
          lm.greens.conservative, lm.greens.progressive)

green.coef <- c(lm.greens.progressive.clustered[2,1], lm.greens.conservative.clustered[2,1],
                lm.greens.market.clustered[2,1],lm.greens.welfare.clustered[2,1])

green.sd <- c(lm.greens.progressive.clustered[2,2], lm.greens.conservative.clustered[2,2],
              lm.greens.market.clustered[2,2],lm.greens.welfare.clustered[2,2])


##
lm.radical_left.welfare <- lm(radical_left~ welfare_chauvinists_voters + 
                                gender + age + as.factor(income) + high.ed + rel.or + unions  + as.factor(countries), 
                              weight=weights, 
                              data=countries.radical.left)
lm.radical_left.welfare.clustered <- coeftest(lm.radical_left.welfare, 
                                              cluster.vcov(lm.radical_left.welfare, countries.radical.left$countries))

lm.radical_left.market <- lm(radical_left~ market_cosmopolitans_voters  +  
                               gender + age + as.factor(income) + high.ed + rel.or + unions  + as.factor(countries), 
                             weight=weights, 
                             data=countries.radical.left)
lm.radical_left.market.clustered <- coeftest(lm.radical_left.market, 
                                             cluster.vcov(lm.radical_left.market, countries.radical.left$countries))

lm.radical_left.conservative <- lm(radical_left~ conservative_voters  +  
                                     gender + age + as.factor(income) + high.ed + rel.or + unions  + as.factor(countries), 
                                   weight=weights, 
                                   data=countries.radical.left)
lm.radical_left.conservative.clustered <- coeftest(lm.radical_left.conservative, 
                                                   cluster.vcov(lm.radical_left.conservative, countries.radical.left$countries))

lm.radical_left.progressive <- lm(radical_left~ progressive_voters  +  
                                    gender + age + as.factor(income) + high.ed + rel.or + unions  + as.factor(countries), 
                                  weight=weights, 
                                  data=countries.radical.left)
lm.radical_left.progressive.clustered <- coeftest(lm.radical_left.progressive, 
                                                  cluster.vcov(lm.radical_left.progressive, countries.radical.left$countries))

stargazer(lm.radical_left.welfare.clustered, lm.radical_left.market.clustered,
          lm.radical_left.conservative.clustered, lm.radical_left.progressive.clustered,
          omit="countries", no.space=TRUE,
          covariate.labels = c("Welfare chauvinists", "Market cosmopolitans",  "Consistent conservatives", "Consistent progressives",
                               "Gender (1=female)", "Age", "Income: medium", "Income: high",
                               "Higher education (=1)",  "Member in religious organization (=1)", "Union member (=1)"))

stargazer(lm.radical_left.welfare, lm.radical_left.market,
          lm.radical_left.conservative, lm.radical_left.progressive)

rl.coef <- c(lm.radical_left.progressive.clustered[2,1], lm.radical_left.conservative.clustered[2,1],
             lm.radical_left.market.clustered[2,1],lm.radical_left.welfare.clustered[2,1])

rl.sd <- c(lm.radical_left.progressive.clustered[2,2], lm.radical_left.conservative.clustered[2,2],
           lm.radical_left.market.clustered[2,2],lm.radical_left.welfare.clustered[2,2])

####Figure 3######
text <- c("CP", "CC", "MC", "WC")

pdf(file='rr_coef.pdf') 
par(mar = c(5.1, 2.1, 4.1, 2.1))
plot(NA, xlim = c(-0.2, 0.2), ylim = c(0, 5), xlab = "", ylab = "", yaxt = "n",
     cex.axis=1.5)
abline(v = 0, col = "gray")
text(x=rr.coef, y=c(1:4), labels=text, cex=2)
segments((rr.coef - qnorm(0.975)*rr.sd)[1], 1, (rr.coef + qnorm(0.975)*rr.sd)[1], 1, col = "gray", lwd = 2)
segments((rr.coef - qnorm(0.975)*rr.sd)[2], 2, (rr.coef + qnorm(0.975)*rr.sd)[2], 2, col = "gray", lwd = 2)
segments((rr.coef - qnorm(0.975)*rr.sd)[3], 3, (rr.coef + qnorm(0.975)*rr.sd)[3], 3, col = "gray", lwd = 2)
segments((rr.coef - qnorm(0.975)*rr.sd)[4], 4, (rr.coef + qnorm(0.975)*rr.sd)[4], 4, col = "gray", lwd = 2)
dev.off()

pdf(file='cd_coef.pdf') 
par(mar = c(5.1, 2.1, 4.1, 2.1))
plot(NA, xlim = c(-0.2, 0.2), ylim = c(0, 5), xlab = "", ylab = "", yaxt = "n",      cex.axis=1.5)
abline(v = 0, col = "gray")
text(x=cd.coef, y=c(1:4), labels=text, cex=2)
segments((cd.coef - qnorm(0.975)*cd.sd)[1], 1, (cd.coef + qnorm(0.975)*cd.sd)[1], 1, col = "gray", lwd = 2)
segments((cd.coef - qnorm(0.975)*cd.sd)[2], 2, (cd.coef + qnorm(0.975)*cd.sd)[2], 2, col = "gray", lwd = 2)
segments((cd.coef - qnorm(0.975)*cd.sd)[3], 3, (cd.coef + qnorm(0.975)*cd.sd)[3], 3, col = "gray", lwd = 2)
segments((cd.coef - qnorm(0.975)*cd.sd)[4], 4, (cd.coef + qnorm(0.975)*cd.sd)[4], 4, col = "gray", lwd = 2)
dev.off()

pdf(file='cons_coef.pdf') 
par(mar = c(5.1, 2.1, 4.1, 2.1))
plot(NA, xlim = c(-0.2, 0.2), ylim = c(0, 5), xlab = "", ylab = "", yaxt = "n",      cex.axis=1.5)
abline(v = 0, col = "gray")
text(x=cons.coef, y=c(1:4), labels=text, cex=2)
segments((cons.coef - qnorm(0.975)*cons.sd)[1], 1, (cons.coef + qnorm(0.975)*cons.sd)[1], 1, col = "gray", lwd = 2)
segments((cons.coef - qnorm(0.975)*cons.sd)[2], 2, (cons.coef + qnorm(0.975)*cons.sd)[2], 2, col = "gray", lwd = 2)
segments((cons.coef - qnorm(0.975)*cons.sd)[3], 3, (cons.coef + qnorm(0.975)*cons.sd)[3], 3, col = "gray", lwd = 2)
segments((cons.coef - qnorm(0.975)*cons.sd)[4], 4, (cons.coef + qnorm(0.975)*cons.sd)[4], 4, col = "gray", lwd = 2)
dev.off()

pdf(file='libs_coef.pdf') 
par(mar = c(5.1, 2.1, 4.1, 2.1))
plot(NA, xlim = c(-0.2, 0.2), ylim = c(0, 5), xlab = "", ylab = "", yaxt = "n",      cex.axis=1.5)
abline(v = 0, col = "gray")
text(x=lib.coef, y=c(1:4), labels=text, cex=2)
segments((lib.coef - qnorm(0.975)*lib.sd)[1], 1, (lib.coef + qnorm(0.975)*lib.sd)[1], 1, col = "gray", lwd = 2)
segments((lib.coef - qnorm(0.975)*lib.sd)[2], 2, (lib.coef + qnorm(0.975)*lib.sd)[2], 2, col = "gray", lwd = 2)
segments((lib.coef - qnorm(0.975)*lib.sd)[3], 3, (lib.coef + qnorm(0.975)*lib.sd)[3], 3, col = "gray", lwd = 2)
segments((lib.coef - qnorm(0.975)*lib.sd)[4], 4, (lib.coef + qnorm(0.975)*lib.sd)[4], 4, col = "gray", lwd = 2)
dev.off()

######Figure A1#####
pdf(file='rl_coef.pdf') 
plot(NA, xlim = c(-0.2, 0.2), ylim = c(0, 5), xlab = "", ylab = "", yaxt = "n",      cex.axis=1.5)
abline(v = 0, col = "gray")
text(x=rl.coef, y=c(1:4), labels=text, cex=2)
segments((rl.coef - qnorm(0.975)*rl.sd)[1], 1, (rl.coef + qnorm(0.975)*rl.sd)[1], 1, col = "gray", lwd = 2)
segments((rl.coef - qnorm(0.975)*rl.sd)[2], 2, (rl.coef + qnorm(0.975)*rl.sd)[2], 2, col = "gray", lwd = 2)
segments((rl.coef - qnorm(0.975)*rl.sd)[3], 3, (rl.coef + qnorm(0.975)*rl.sd)[3], 3, col = "gray", lwd = 2)
segments((rl.coef - qnorm(0.975)*rl.sd)[4], 4, (rl.coef + qnorm(0.975)*rl.sd)[4], 4, col = "gray", lwd = 2)
dev.off()

pdf(file='green_coef.pdf') 
plot(NA, xlim = c(-0.2, 0.2), ylim = c(0, 5), xlab = "", ylab = "", yaxt = "n",      cex.axis=1.5)
abline(v = 0, col = "gray")
text(x=green.coef, y=c(1:4), labels=text, cex=2)
segments((green.coef - qnorm(0.975)*green.sd)[1], 1, (green.coef + qnorm(0.975)*green.sd)[1], 1, col = "gray", lwd = 2)
segments((green.coef - qnorm(0.975)*green.sd)[2], 2, (green.coef + qnorm(0.975)*green.sd)[2], 2, col = "gray", lwd = 2)
segments((green.coef - qnorm(0.975)*green.sd)[3], 3, (green.coef + qnorm(0.975)*green.sd)[3], 3, col = "gray", lwd = 2)
segments((green.coef - qnorm(0.975)*green.sd)[4], 4, (green.coef + qnorm(0.975)*green.sd)[4], 4, col = "gray", lwd = 2)
dev.off()

pdf(file='socdem_coef.pdf') 
plot(NA, xlim = c(-0.2, 0.2), ylim = c(0, 5), xlab = "", ylab = "", yaxt = "n",      cex.axis=1.5)
abline(v = 0, col = "gray")
text(x=socdem.coef, y=c(1:4), labels=text, cex=2)
segments((socdem.coef - qnorm(0.975)*socdem.sd)[1], 1, (socdem.coef + qnorm(0.975)*socdem.sd)[1], 1, col = "gray", lwd = 2)
segments((socdem.coef - qnorm(0.975)*socdem.sd)[2], 2, (socdem.coef + qnorm(0.975)*socdem.sd)[2], 2, col = "gray", lwd = 2)
segments((socdem.coef - qnorm(0.975)*socdem.sd)[3], 3, (socdem.coef + qnorm(0.975)*socdem.sd)[3], 3, col = "gray", lwd = 2)
segments((socdem.coef - qnorm(0.975)*socdem.sd)[4], 4, (socdem.coef + qnorm(0.975)*socdem.sd)[4], 4, col = "gray", lwd = 2)
dev.off()

########table A12#######
uk.1990 <- data.1990[data.1990$countries=="Great Britain",]

uk.right_b.delta.90 <- lm(binary.right~red.minus.values + I(red.minus.values^2), 
                          weight=weights,
                          data=uk.1990)
summary(uk.right_b.delta.90)

uk.right_b.delta.covariates.90 <- lm(binary.right~red.minus.values + I(red.minus.values^2) + 
                                       gender + age + as.factor(income) + rel.or + unions, 
                                     weight=weights,
                                     data=uk.1990)
summary(uk.right_b.delta.covariates.90)

uk.2017 <- data.2017[data.2017$countries=="Great Britain",]

uk.right_b.delta.17 <- lm(binary.right~red.minus.values + I(red.minus.values^2), 
                          weight=weights,
                          data=uk.2017)
summary(uk.right_b.delta.17)
uk.right_b.delta.covariates.17 <- lm(binary.right~red.minus.values + I(red.minus.values^2) + 
                                       gender + age + as.factor(income) + high.ed + rel.or + unions, 
                                     weight=weights,
                                     data=uk.2017)
summary(uk.right_b.delta.covariates.17)

netherlands.1990 <- data.1990[data.1990$countries=="Netherlands",]

netherlands.right_b.delta.90 <- lm(binary.right~red.minus.values + I(red.minus.values^2), 
                                   weight=weights,
                                   data=netherlands.1990)
summary(netherlands.right_b.delta.90)

netherlands.right_b.delta.covariates.90 <- lm(binary.right~red.minus.values + I(red.minus.values^2) + 
                                                gender + age + as.factor(income) + rel.or + unions, 
                                              weight=weights,
                                              data=netherlands.1990)
summary(netherlands.right_b.delta.covariates.90)

netherlands.2017 <- data.2017[data.2017$countries=="Netherlands",]

netherlands.right_b.delta.17 <- lm(binary.right~red.minus.values + I(red.minus.values^2), 
                                   weight=weights,
                                   data=netherlands.2017)
summary(netherlands.right_b.delta.17)

netherlands.right_b.delta.covariates.17 <- lm(binary.right~red.minus.values + I(red.minus.values^2) + 
                                                gender + age + as.factor(income) + high.ed + rel.or + unions, 
                                              weight=weights,
                                              data=netherlands.2017)
summary(netherlands.right_b.delta.covariates.17)

stargazer(uk.right_b.delta.90, uk.right_b.delta.covariates.90,
          uk.right_b.delta.17, uk.right_b.delta.covariates.17, 
          netherlands.right_b.delta.90, netherlands.right_b.delta.covariates.90,
          netherlands.right_b.delta.17, netherlands.right_b.delta.covariates.17,
          no.space=TRUE,
          covariate.labels = c("Delta EC", "Delta EC^2", "Gender (1=female)", "Age", "Income: medium", "Income: High",
                               "Higher education (=1)",  "Member in religious organization (=1)", "Union member (=1)"))

######Table A13####
right_b.delta.interest.90 <- lm(binary.right~red.minus.values + I(red.minus.values^2) + 
                                  gender + age + as.factor(income) + rel.or + unions + as.factor(interest) + as.factor(countries), 
                                weight=weights,
                                data=data.1990)
right_b.delta.interest.90.clustered <- coeftest(right_b.delta.interest.90, 
                                                cluster.vcov(right_b.delta.interest.90, data.1990$countries))

right_b.delta.interest.17 <- lm(binary.right~red.minus.values + I(red.minus.values^2) + 
                                  gender + age + as.factor(income) + high.ed + as.factor(interest)+ rel.or + unions + as.factor(countries), 
                                weight=weights, 
                                data=data.2017)
right_b.delta.interest.17.clustered <- coeftest(right_b.delta.interest.17, 
                                                cluster.vcov(right_b.delta.interest.17, data.2017$countries))

stargazer(right_b.delta.interest.90.clustered,right_b.delta.interest.17.clustered,
          omit="countries", no.space=TRUE,
          covariate.labels = c("Delta EC", "Delta EC^2", "Gender (1=female)", "Age", "Income: medium", "Income: High",
                               "Higher education (=1)",  "Member in religious organization (=1)", "Union member (=1)",
                               "Interest: somewhat interested", "Interest: not very interested", "Interest: not at all interested"))

stargazer(right_b.delta.interest.90, right_b.delta.interest.17)

#####Table A14#####
right_b.delta.truncated.90 <- lm(binary.right~red.minus.values.truncated + I(red.minus.values.truncated^2) + 
                                   gender + age + as.factor(income) + rel.or + unions + as.factor(countries), 
                                 weight=weights,
                                 data=data.1990)
right_b.delta.truncated.90.clustered <- coeftest(right_b.delta.truncated.90, 
                                                 cluster.vcov(right_b.delta.truncated.90, data.1990$countries))

right_b.delta.truncated.17 <- lm(binary.right~red.minus.values.truncated + I(red.minus.values.truncated^2) + 
                                   gender + age + as.factor(income) + high.ed + rel.or + unions  + as.factor(countries), 
                                 weight=weights, 
                                 data=data.2017)
right_b.delta.truncated.17.clustered <- coeftest(right_b.delta.truncated.17, 
                                                 cluster.vcov(right_b.delta.truncated.17, data.2017$countries))

stargazer(right_b.delta.truncated.90.clustered, right_b.delta.truncated.17.clustered,
          omit="countries", no.space=TRUE,
          covariate.labels = c("Delta EC", "Delta EC^2", "Gender (1=female)", "Age", "Income: medium", "Income: High",
                               "Higher education (=1)",  "Member in religious organization (=1)", "Union member (=1)"))
stargazer(right_b.delta.truncated.90, right_b.delta.truncated.17)

#####Table A15######

right_no_middle.delta.covariates.90 <- lm(lr_no_middle~red.minus.values + I(red.minus.values^2) + 
                                            gender + age + as.factor(income) + rel.or + unions  + as.factor(countries), 
                                          weight=weights, 
                                          data=data.1990)
right_no_middle.delta.covariates.90.clustered <- coeftest(right_no_middle.delta.covariates.90, 
                                                          cluster.vcov(right_no_middle.delta.covariates.90, data.1990$countries))

right_no_middle.delta.covariates.17 <- lm(lr_no_middle~red.minus.values + I(red.minus.values^2) + 
                                            gender + age + as.factor(income) + high.ed + rel.or + unions  + as.factor(countries), 
                                          weight=weights, 
                                          data=data.2017)
right_no_middle.delta.covariates.17.clustered <- coeftest(right_no_middle.delta.covariates.17, 
                                                          cluster.vcov(right_no_middle.delta.covariates.17, data.2017$countries))

stargazer(right_no_middle.delta.covariates.90.clustered,right_no_middle.delta.covariates.17.clustered,
          omit="countries", no.space=TRUE,
          covariate.labels = c("Delta EC", "Delta EC^2", "Gender (1=female)", "Age", "Income: medium", "Income: High",
                               "Higher education (=1)",  "Member in religious organization (=1)", "Union member (=1)"))
stargazer(right_no_middle.delta.covariates.90, right_no_middle.delta.covariates.17)

#########Table A1#######
rm(list = ls())

load(file='~/Dropbox/many ways to be right submission/BJPS final submission/BJPS replication materials/data_2017_factor.RData')
load(file='~/Dropbox/many ways to be right submission/BJPS final submission/BJPS replication materials/data_1990_factor.RData')

factor.1990 <- fa(recoded.1990.factor, nfactors=2, rotate="promax",scores=T)
loading.1990 <- factor.1990$loadings
loading.1990[,c(1,2)]
#SS loadings==1.06, 1

factor.2017 <- fa(recoded.2017.factor, nfactors=2, rotate="promax",scores=T)
loading.2017 <- factor.2017$loadings
loading.2017[,c(1,2)]
#SS loadings==1.33, 1.1

ss.loadings <- c(1.06, 1, 1.33, 1.01)

table_factor <- cbind(loading.1990[,c(1,2)],loading.2017[,c(1,2)])
table_factor <- rbind(table_factor, ss.loadings)
xtable(table_factor)
